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ABSTRACT 



The origin of hard X-ray (HXR) excess emission from clusters of galaxies is 
still an enigma, whose nature is debated. One of the possible mechanism to pro- 
duce this emission is the bremsstrahlung model. However, previous analytical 
and numerical calculations showed that in this case the intracluster plasma had 
to be overheated very fast because suprathermal electrons emitting the HXR ex- 
cess lose their energy mainly by Coulomb losses, i.e., they heat the background 
plasma. It was concluded also from these investigations that it is problematic 
to produce emitting electrons from a background plasma by stochastic (Fermi) 
acceleration because the energy supplied by external sources in the form of Fermi 
acceleration is quickly absorbed by the background plasma. In other words the 
Fermi acceleration is ineffective for particle acceleration. We revisited this prob- 
lem and found that at some parameter of acceleration the rate of plasma heating 
is rather low and the acceleration tails of non-thermal particles can be gener- 
ated and exist for a long time while the plasma temperature is almost constant. 
We showed also that for some regime of acceleration the plasma cools down in- 
stead of being heated up, even though external sources (in the form of external 
acceleration) supply energy to the system. The reason is that the acceleration 
withdraws effectively high energy particles from the thermal pool (analogue of 
Maxwell demon). 

Subject headings: galaxies: clusters: individual (Coma) — X-rays — physical data 
and processes 
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Introduction 



One of the most important problem in astrophysics is the problem of particle 
acceleration. The general expression for acceleration of a charged particle is 



y(7Wv) = Ze ( E + -v x B 
at \ c 



(1) 



where E and B are the electric and magnetic field strength, v is the velocity of particle 
and 7 = 1/ a/1 — v 2 /c 2 . In most astrophysical conditions static electrical fields cannot be 
maintained because of a very high electrical conductivity. Therefore the acceleration can 
be associated either with non-stationary electric fields (electromagnetic waves) or with 
time-varying magnetic fields. In the latter case the work can be done by the induced electric 
field 

-VxE. (2) 



ldB 

c~dt 



The basic idea of accelerati on by 



conditions was suggested by 



electromagn etic inhomogeneities in astrophysical 



Fermil (11949 



1954 ) who assumed that the Galactic cosmic rays 



(CRs) were accelerated by collisions of charged particles with fluctuations of magnetic fields 
(magnetic clouds) moving chaotically with the velocity dispersion u. One of the features of 
this theory was that it yielded naturally a power-law spectrum of accelerated particles. The 
rate of particle acceleration by this stochastic mechanism is about 



\dt) F 



u f 
~ — I 



(3) 



where £ and v are the particle kinetic energy and velocity, and r is the average time of 
particle collision with the clouds. This rate of acceleration is slow because u<C«. In kinetic 
equations, the Fermi (stochastic) acceleration is described as momentum diffusion (see e.g., 



Toptygin 



19851) 



df_ _ }_d_ 

dt p 2 dp 



D F (p)p 



2 df 
dp 



+ Lf = 



(4) 
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with the diffusion coefficient Dp(p) in the form 



D F (p) ~ p" 



u 

V 2 T 



(5) 



Here f(p,t) is the particle distribution function, p is the particle momentum, t is the time, 
and the operator L describes particle spatial propagation and their momentum losses. 

In spite of its low efficiency, st ochastic accelerat i on may be esse ntial for particle 



acceleration in solar flar es (see e.g., 



medium of the Galaxy 



Miller et al. 



Mertsch fc Sarkar 



201l|). 



( jBerezinskii et al 



1990; 



Petrosian 



20121 ). in the interstellar 



19901 ) and near the Galactic center (see e.g., 



The problem of stochastic particle acceleration in galax y clusters arose from 

HXR) energy range (see e.g., 



observations in the hard X-ray 



2007 



Rephaeli et al. 



1999 



2008 



Eckert et al. 



2008 



Fusco-Femiano et a 



Nevalainen 



2009 



Ajello et al. 



1999, 



20101 ) 



which showed an emission excess above the equilibrium thermal X-ray spectrum. 

One of the several interpretations of the HXR excess from the Coma cluster in the 
range 20-80 keV was an assump tion that it was pr oduced by bremsstrahlung radiation of 



suprathermal electrons (see e.g., 



However, this model was criticized by 



Petrosian 



Enfilin et al.lll999f) accelerated in the intracluster medium. 



(|200ll ) who concluded from simple estimates 



that in this case the intracluster plasma in Coma had to be overheated very fast. The 
point is that suprathermal electrons emitting the HXR excess lose their energy mainly by 
Coulomb losses, i.e., they lose their energy by heating the background plasma. If these 
electrons generate an X-ray flux Lx by bremsstrahlung, they transfer the energy flux Lc to 
the background plasma. The necessary energy input is estimated as 

(dS/dt) c 



■>c 



(d£/dt) B R 

where (dE/dt)c arid (dS/dt)sR are the rates of Coulomb and bremsstrahlung losses, 
respectively In the keV energy range (d£/dt) c (d£/dt) BR , and this seems to make 



(6) 
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plasma overheating inevitable. However, we should point out that in effect particle 
acceleration may also be accompanied by plasma cooling due to run-away flux of high 
energy particles from thermal pool, and more careful analysis is necessary to define which of 
these effects (plasma heating or cooling) prevails. This analysis is presented in the following 
sections. 



2. Review of Particle Acceleration from Background Plasma 

A natural source for suprathermal particles is stochastic acceleration of seed particles 
from a background plasma. These particles are accelerated when the rate of acceleration 
(dE/d£)p exceeds the rate of the Coulomb losses (d£/dt)c- A characteristic energy called 
the injection energy £ m j is the energy above which a non-thermal spectrum is formed by 
acceleration. It is determined by equating these rates of acceleration and loss. 

The kinetic equation in the particle momentum space describing stochastic particle 
acceleration from background plasma has the form (assume isotropic distribution) 



dt p 2 dp 



0, (7) 



where Dp(p) is the diffusion coefficient of stochastic (Fermi) acceleration, and (dp/dt)c and 
Dq{p) describe particle momentum losses and diffusion due to Coulomb collisions. These 
coefficients are calculated from the total distribution function / (see Appendix |XJ) , and 
therefore in general the equation is nonlinear. 



2.1. Linear Approximation 



An analytical solution of this equation for the case of weak a c celera tion from a 
background plasma with temperature T was obtained by iGurevichl (119601 ) . The term of 



-6 - 



stochastic acceleration was taken in the phenomenological form 



D F {p) = ap 2 . (8) 

The analysis was provided for the case when the characteristic time of stochastic 
acceleration, 

T F =p 2 /D F} (9) 

is much larger than the time of thermal particle collisions, r t h, 

~2 m e {k B Tfl 2 



Tth 



(10) 



m vriVe 4 lnA ' 

where N is the density and T is the temperature of background plasma, In A is the Coulomb 
logarithm, m e is the electron rest mass and m is the mass of accelerated particles. 

In this case the injection energy £ in j is much larger than the plasma temperature, 
£i n \ ^> k B T. Coulomb collisions keep the equilibrium Maxwellian distribution for most 



part of the momentum ran ge, and t 



p ^> \/2mk B T (as used by 



re coe fficients of Eq. ([7]) for nonrelativistic momenta 



Gurevich 



19601 ) for the Maxwellian distribution function. For 



T th t f , significant distortions from the equilibrium Maxwellian state are expected only for 
very large values of mome nta and a very small fraction of thermal particles is accelerated. 
Therefore 



Gurevich! (119601 ) assumed that the number of particles N(t) in the momentum 
range p < p u ^ varies very slowly with time t, N(t) — Nq — St, where Nq is the initial particle 
density and a small run-away flux S is generated at relatively high momentum range. The 

run-away flux for the case of slow acceleration can be described as 

P 



S{p) = S c 



x 2 e x dx = Sn 



erf (p) — p e 



(11) 



where erf(z) is the error function, p = pj y/2mk&T, and the constant So is derived from 
boundary conditions. The flux is zero at p = but when p 3> \/2mk B T, it reaches a 
maximum value S(p) = So as shown in Fig. [TJ 
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Fig. 1. — Run-away flux S as a function of the dimensionless momentum pj \/1mk^T . 



In the momentum range where S(p) ~ 5*0 the distribution function is non-Maxwellian 
and it is described by the kinetic equation 



|)^_ {DcW+flp(p)} | 



c 



S . 



(12) 



The acceleration forms a nonthermal component of the spectrum in the range p > p iri j where 
p in j is the solution of equation 

However, if D F (p) ^ in the range p < p in j, a solution of this equation describes also 
an excess of the distribution function above the equilibrium Ma xwellian distribution in 



Gurevich 



1960). This excess at p < p in j is 



momentum ranges both above and below p in j (see 
formed by Coulomb collisions in the transition range between the thermal (Maxwellian) and 
non-thermal parts of the spectrum. If the HXR excess is due to b remsstrahlung em ission of 



electrons from this transition region then the relation ([6]) used by 



Petrosian 



(120011 ) cannot 



be applied to the estimate of Lc and more accurate calculations are necessary. 



Br emsst r ahlung emiss i on of e 



Dogiell fl2000h : 



Liang et al 



J2002I); 



ectrons from the transition region was calculated in 

The conclusion is that the necessary 



Dogiel et al. 



(120071 ). 
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energy input Ln for Coma was about one order of magnitude less than obtained by 



Petrosianl ( 1200 ll ). This may solve the problem of the plasma overheating. However, their 



linear analysis of equation (J7J) does not include variations of temperature T which is 
supposed to be constant. 



2.2. Non-Linear Treatment 



More reliable conclusions ca n be derived from analys es of the nonlinear equation in 



the form similar to those used by 



MacDonald et al. 



(jl957l ). when a feedba ck of accelerate d 



partic les o n the plasma temp e ratur e is taken into account. Very recently 



Wolfe fc Melia 



( j200fih and 



Petrosian East 



( 120081 ) provided similar numerical analysis for the case of 



stochastic acceleration from a background plasma. 

The nonlinear kineti c equation describing particle Coulo mb collisions is deriv e d in 



Landau fc Lifshitzl ( I198ll ) (see Appendix [SJ. Using this theory 



Navakshin fc Melial (119981 ) 



derived coefficients of this equation for the case of isotropic and h omogeneous distributi on 



function for non-relativistic and ultra-relativistic particles. Later, 



Wolfe k Melial (120061 ) 



extended their analysis to the general case of anisotropic distribution function. 



Numerical analysis of these equations has been performed by 
the isotropic stochastic acceleration in the form 

D F (p) = ap"6(p - 1/2) . 



Wolfe fc Melial (|2006|) for 



(14) 



Wolfe fc Melial (120061 ) stated that the continuous stochastic acceleration of thermal electrons 
produced a nonthermal tail. But for the hard X-ray emission in the Coma Cluster this 
model actually cannot work because the energy gained by the particles is distributed to 
the whole plasma on a timescale much shorter than that of the acceleration process itself. 
Moreover, bremsstrahlung is relatively inefficient to cool the accelerated electrons, the 
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energy of this tail is quickly dumped into the thermal background plasma and heat the 
plasma. 



Similarly, iPetrosian fc Eastl (120081 ) obtained numerical solutions of the nonlinear 
isotropic kinetic equations which included effects of plasma heating for the stochastic 
diffusion in the form 



D F (S) 



S 2 



((£)t (1 + £ c /£)« 



(15) 



where £ = ^/p 2 + 1 — 1 is the kinetic energy normalized to mc 2 , C(£) = (2 — 7 2 )/(l + 7 1 ), 
and To, £ c and q are free parameters. 



Petrosian fc East 



results of 



Dogiel et al 



(2008) concluded that their calculations confirmed qualitatively 
(120071 ) that the required input energy Lc was lower than that follows 
from the estimate but by a factor of 2 or 3 only, and that did not solve the problem 
of plasma over h eating . Besides, they argued that their calculations confirmed results of 



Wolfe fc Melial (120061 ) that stochastic acceleration could not work in clusters because the 
energy gained by the particles was distributed to the whole plasma on timescales much 
shorter than that of the acceleration process. At acceleration rates smaller than the 
thermalization rate of the background plasma, there is very little acceleration. The primary 
effect of acceleration is heating of the plasma. In the opposite case, at higher energizing 
rates, a distinguishable nonthermal tail is developed, but this is again accompanied by an 
unacceptably high rate of heating. 

In other words it follows from these investigations that it is problematic to accelerate 
particles from a background plasma because the main effect of this acceleration is plasma 
overheating. The energy supplied by external sources in the form of stochastic (Fermi) 
acceleration is quickly absorbed by a background plasma. An interesting question 
arises: whether any conditions exist when the stochastic acceleration generate prominent 
nonthermal tails while the plasma is not overheated and its temperature varies relatively 
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slowly. From the analysis in the following sections, we argue that the answer is affirmative. 



3. Particle Acceleration from Background Plasma: Quasi-Linear 

Approximation 

First, we estimate variations of plasma temperature derived in quasi-stationary 
approximations when the distribution function can be presented as f — f(p,N,T). In this 
case, 

df = df_dN d£dT 

dt ON dt dT dt ' 1 ' 

where N = N(t) and T = T(t) are slowly varying functions of t. 



3.1. Distribution function 

In this subsection we investigate the isotropic form of the kinetic equation (lAlj) . This 
equation describes stochastic particle acceleration from background plasma and it is exactly 
the same as Eq. (j7|). The appropriate boundary conditions are Eqs. flA5[) & (IA6j) . Recall 
that the particle momentum has been normalized to mc. Here and in the following the 
temperature T is indeed the thermal energy k^T normalized to mc 2 . The particle kinetic 
energy 8 = ^/p 2 + 1 — 1 is also normalized to mc 2 . The coefficients {dp/dt)c{p), Dc{p) and 
Dp(p) are normalized accordingly. 

The stochastic Fermi acceleration is supposed to be isotropic and has a phenomenological 
form as 

D F (p) = ap"6(p-p ), (17) 
where a, q and po are arbitrary parameters. The problem is characterized also by the 
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injection momentum 



tttfnj = -Pinj ( ^ 



;is) 



P=Pinj 

The acceleration is effective in the momentum range p > max{p ,Pinj}- 



Similar to iGurevichl (jl960f ) we assume that the acceleration time tf, is much longer 
than the time of thermal particle collisions r t h, i.e., values of p in j or po are large and one of 
the corresponding energy values is much higher than the temperature, 

T<max(£ inj ,£ ). (19) 

In this case, Coulomb collisions keep the equilibrium Maxwellian distribution over an 
extended momentum range with a significant deviation from this distribution at very large 
momenta, i.e., a small part of thermal particles is accelerated. The number of non-thermal 
particles generated by the acceleration N n in this case is much smaller than the number of 
thermal particles N, N n /N <C 1. 

Below we present the distribution function and the coefficients of the kinetic equation 
as series expansions over the small parameter e = N n /N <C 1, 

f(p,t) = fo(p,t) + h(p,t) + 0(e 2 ) , 
D c (p,t) = D (p,t) + D 1 (p,t) + O(e 2 ) , (20) 



Here O (e l ) denotes terms of order e l or above. Note that fi(p,t) = 0(e l ), D and (dp/dt)o 
are calculated from Eq. (1A2j) for the function / , and D x and (dp/dt)i for the function / 1( 
etc. 

In the quasi-stationary approximation the derivative df/dt can be presented in the 
form ( [16]) . The derivatives dN/dt and dT/dt can be presented as series dN/dt = 0(e) and 
dT/dt = 0(e), because without acceleration (iV n = 0) we have dN/dt = and dT/dt = 0. 
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Here we have 



and dfo/dt is of the order of e. 



d£ = dfo 
dt dt 



+ 0(e 2 ) 



(21) 



It is convenient to express the distribution function as 



f(p) = f(p)8(po -p) + f H (p)e(p - pa) 



(22) 



First, we find the solution of Eq. ([7]) in the momentum range < p < po where the 
acceleration term vanishes and / = / (see Eq. (1221) ). In zero order of expansion (no 
acceleration) the function /o is Maxwellian 

rp / dp\ dp 



fo(p) = Coexp 



dt J D 



C exp(-£/T) 



(23) 



where D and (dp/dt)o are the Maxwellian kinetic coefficients. For p 3> a/T 2 + 1 — 1 the 
Bethe-Bloch approximation for the these coefficients is 

1 



dp 
lit 



-A 1 + 



P 



(24) 



-Tjl 



1 f dp 



p 2 \dt 



AT 1 



1 

P 2 



3/2 



Here and below 



A = 4nKcN In A . 



(25) 



The characteristic time of Coulomb losses for a particle with momentum p (in unit of 
mc) is 



Tcip) 



P 



A(p 2 + 1) 



(26) 



Constant Cn is estimated from the normalization condition 



C = N 



Po 



P 2 fo(p) d P 



-1 



iVexpt-T- 1 ; 
TK 2 {T- 1 ) 



(27) 
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where K2(x) is the modified Bessel function. For non-relativistic temperatures (T <C 1) we 
obtain 

(2? 



The kinetic equation for the function // can be rewritten as 



1 9 2 

^2 



9p 



<9p 



dt 



at 



0(e 2 ) 



(29) 



Integrating the above equation from to p gives 



P 



D (p) 



df[ 
dp 



dp) f i 
dtJ J \ 



S = - (S 1 + S 2 ) 



(30) 



where S is the flux of particles through the point p. Here 



Si 



dN(p, t) 
dt 



S 2 = P 2 



p 

d f 

■qI / u 2 fv(u)du, 



(31) 



(32) 



The flux Si describes a particle leakage (in momentum space) caused by the acceleration. 
It generates a slow decrease of particle number in the thermal region. The flux S 2 causes 
the plasma heating and temperature variations with time. 



Thus, the solution of Eq. (I3"U1) is 



flip) = exp(-£/T) 



As in 



Ci 



S(u) 



u 2 D (u) 



exp(£ /T)du 



(33) 



Gurevichl ( 119601 ) the value of the constant C\ can be derived from the normalization 
condition 

PO 

Jp 2 fi(p)dp = 0. (34) 


The kinetic coefficients D\ and {dpjd£)\ are calculated for the function f( + f 11 . Therefore 
Eq. f )33|) is an integral equation for f{{p) which should be added by an equation for f H {p). 
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The asymptotic form of flip) for large values of p can easily be derived. Indeed, if £ ^> T 
then Slip) = 0(e) while S 2 ip) ~ 0(e) exp(— £ /T) <C S\. Therefore S 2 ip) can be neglected. 
As one can see from Eq. f[3~Tj) the flux Si(p) remains almost constant for sufficiently large p 
(see Fig. [p. So w ith a high degree of accuracy we can put Slip) = Sn = —dN/dt, which is 



the same as So in 



Gurevichl (Jl960|) 



It follows from Eqs. f )33|) and ( l3~4j) that the constant C\ is 

O 1 «^ JV T C (p )y|T- 3 / 2 , 



(35) 



where Tcipo) is the characteristic time of Coulomb collision for the particle momentum 
p = Po (for t f and r c see Eqs. OH]) and (1261) ). Thus, for the estimation of obtained in 
subsection 13.31 we have 



iV f £ \ r c jp ) N ( £ 
The distribution function in the range p < po can be written as 



(36) 



fiv) * tiip) + flip) 
N 



TKziT- 1 ) 6XP V T7 ,17 



5 



A 1 



(37) 



where i = \f p 2 + 1 = £ + 1 is the total energy of particle and 

exp(x) 



Eiiz) 



-dx . 



x 



(31 



For non-relativistic temperatures £/T ^> 1 the expansion of Eiiz) for z ^> 1 is 



exp(z 



k=0 



(39) 



Thus for large values of p 

fiv) 



2 N ( £ 

7TT^ eXP 



S 



N 



TJ Aip 2 + V 



(40) 
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The distributions function Eq. ( I37|) can be presented in the form 
f\p) = 



fi(p) + 0(e), ioiS<T 



(41) 



In the range p > po the acceleration cannot be neglected. With the constant flux Sn of 
particles the equation for the distribution function f 11 in this region reads 

df n 



p~ 



{D (p) + D F (p)} 



dp 



The general solution of this equation is (see, e.g. 



fdp\ f ii 
\dt)/ _ 




Gurevich 


1960) 



-s 



N 



(42) 



| P 

f n (p) = C n e W U 



(dp/dt) (u)du 
D F (u) + D (u) 

p 



-S N exp 



v 2 dv 



(dp/dt) (u)du 
D F (u) + D (u) \ J D F (v) + D (v 



■ exp 



(dp/dt) (u)du 
D F {u) + D (u) 



.(43) 



The constant C H can be estimated from the continuity condition at p = p : f T (po) = f n (po) 
while the value of Sn can be estimated from the second boundary condition: / /7 (p m ax) = 0. 

For p 3> Pmj, we can assume acceleration dominates Coulomb loss. It is easy to show 
from Eq. (|42]) and Eq. (fTTJ) that the function f n {p) is a power-law 

Sn 



ell 



(p) =C X + 



p-"- 1 , 



(44) 



where G\ is a constant. 



3.2. Plasma heating rate 

Using the total distribution function / (see Eq. (1221) . where f 1 and f 11 are determined 
by Eqs. (|4T|) & f|43l ). we can calculate the kinetic coefficients (1A2I) for the nonlinear 
equation ([7]), and then estimate the temperature variations of the background plasma 
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caused by particle acceleration. In this case the stochastic Fermi momentum diffusion 
describes the energy supply into the system by external sources. Generally speaking, energy 
supply can vary with time, but usually it is assumed that external sources keep a stationary 
level of acceleration such that Dp is constant. 



The total energy input into the system is 

oo 

Wext = " I 

J dp 



' 2n d f 
P Dr- 
op 



dp . 



(45) 



It is a function of time even if Dp is constant, because the distributions function / is time 
dependent. 

Note that Coulomb collisions do not change the total energy in the system, therefore 
we have 

oo 

r ft _ r / An \ ft f 

dp = 0. (46) 



S 7rP 
op 



dt J r , dp 



c 



This condition is valid for any function / if the kinetic coefficients (dp/dt)c and Dc(p) are 
calculated from Eq. f ]A2j) for this function /. 



The energy supplied by the stochastic Fermi acceleration is distributed over the 
spectrum in the form of accelerated particles and a heated plasma, because accelerated 
particles lose their energy by Coulomb collisions and thus transfer a part of their energy to 
thermal particles. Variations of dT/dt in the quasi-equilibrium part of the spectrum can be 
derived from estimates of the energy flux into the region p < p which is 



PO PO 

W = -^J P 2 £f 1 (p)dp = j 8 



d_ 

dp 



P 2 D, 



dp 



p 



dp 
~dt 



/ c 



dp. 



(47) 



o o 

where the coefficients D c and (dp/dt) c are calculated for the total distribution function (T52 
For the estimation of the integral f )4T|) we can use the condition (JM]), and obtain 

oo 



Wn 



dp 



P 2 D C - 



dp 



P 



2 I dP \ rll 
dt ' 



dp . 



(48) 



PO 



-re- 



integration by parts gives 

Wo = S S N + 



P 



PO 



dp . 



(49) 



Since Wq = 0(e) and f 11 = 0(e) we can use the Maxwellian (Bethe-Bloch) expressions for 
the kinetic coefficients D and (dp/dt) as in Eq. (p4|) . 

OO 

5/" 



W„ = -£ Sn + 



V 



pit 



dp 



dp 



f 



dp . 



(50) 



We see that the energy input into the thermal part of the spectrum (plasma heating) 
is determined by two processes: (i) energy losses of nonthermal particles (the integral of 
Eq. ( 15 Op ) which heat the plasma, and (ii) a particle escape to the high energy part of the 
equilibrium spectrum (the first term on the RHS of Eq. (15 Op ) which cools the plasma. On 
the other hand we can express Wo in the form 



dW dW dT 8W dN 



dt 



dT dt dN dt ' 



(51) 



To the first order of e of the expansion of dT /dt we can take Wq as 



PO 

Wq = J u 2 £f^(u) du 



PO 



iVexp(-T- 1 ) f 2 / £\ 
——— rr — u t exp — — du . 
TK 2 (T- 1 ) J P V T' 





(52) 



In the general case the temperature variations can be calculated numerically (see section H]). 
However, these calculations can be simplified. The point is that the particle spectrum 
described by Eqs. (HOI) & (jl5]l depends strongly on the relation between the momenta pj n j 
and po. Fig. |2] illustrates this situation: as po increases the transition region in momentum 
range p > po shrinks and finally disappears when po reaches p- m j. In the limiting case 
Po > Pmj the transition region vanishes almost completely and the power-law tail of 
nonthermal particles is attached almost directly to the thermal equilibrium distribution. 
In this case, evaluations of the plasma temperature can be performed analytically because 
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the fu nctions f 1 and f 11 have very simple form. We notice that the conclusion of iGurevich 



( Il960l ) about a very extended transition region between thermal and nonthermal parts of 



the spectrum is valid only for the case when p < p in j. 




Fig. 2. — Shape of f(p) for different values of pq. Thin solid line represents f(p), thick solid 
line - pure Maxwellian distribution, dashed line - power-law approximation. 



3.3. The case of transitionless acceleration 

If Pq > pi n j, Eq. ( jH]) is an appropriate solution for the distribution function. C\ 
and Sn are determined from the boundary conditions at p = po and p = p max , namely, 
f H (po) = f(po) = fo and f n ( Pmax ) = 0, 

S N = a(q + l)p ( +1 fo , (53) 
q „-(?+!) /„ \ -(«+!) 

a{q + 1) V Po J 
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As p m ax ^ Po, thus for simplicity we set C\ = 0, and for non-relativistic temperatures 
T 1 from Eq. (HOI we have 

In this case the run-away particle flux toward high energies can be expressed directly from 
Eq. ([MD as 

S N = a(? + l)pl+^JL exp 
For po ^> 1 Eq. (IBlZj) becomes 

oo 

Wo = y p*Sf*(p)dp = N 
o 

or for non-relativistic values of T <C 1 

3 15 

W = -NT + —NT 2 + ... (58 
2 8 



1 + 



(55) 



1 + 



«(<? + l)Po +1 



(56) 



(3T - 1) + 



K X {T 



K 2 (T-i) 



(57) 



Now we have (recall Eqs. (J50|) and fl5T|) ) 



«9T tZt 
= a/o^oPo +1 (^ + l) 



— - t J &N 



P 



Po 



a£ (? + 1) 



p 2 + 1 L ^ \ tft / o 

f 3«J9 +1 (c; + l) 



V 



(59) 



(^-1) J 



where 



Q(po,?) = y x~Wx 2 + 1 dx . 



(60) 



If a£o{s + 1) 7^ ^4Q(po, ?) and So ^> T then the second term in Eq. fl59|) is small and can be 
neglected. Finally from Eq. ( |58l) we obtain 



v 



dt 3iV 



ak + 1) 



So 



(61) 
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where A and Sn are defined by Eqs. f l25|) and fl56|) . 

For high values of a one can see from Eq. (I6T1) that the plasma cools down and 
dT/dt < 0. The temperature decreases with time due to a very intensive outflow of 
high energy particles from the thermal pool, even though external sources in the form of 
stochastic Fermi acceleration supply energy to the system (analogue to Maxwell demon). 
This effect can be seen in Fig. [2] as a deficit of high energy thermal particles at p < po. 

When a decreases, collisions start to dominate over the outflow effect. As the result 
the derivative dT/dt increases and at sufficiently small a the regime changes from cooling 
to heating of plasma. However the process of acceleration reduces the amount of particles 
in the thermal pool (dN/dt < at p < po). If a remains constant then the value of A 
decreases with time and, in principle for a sufficiently long time we come again to the 
condition when a£o(q + 1) > A(t)Q(p , ^), enter the regime of plasma cooling again. 

A more accurate analysis of this regime can be provided by numerical calculations of 
the nonlinear case. 

4. Nonlinear Case: Semi-Analytical Method and Numerical Calculations 

The most straightforward way to solve the problem is a numerical solution of the 
original nonlinear equation. However, this method is very time-consuming. We proceed 
with approximation methods that simplify the numerical calculations, but still give a good 
result. 

Analysis of kinetic equations depend on the relation between the plasma heating time 
and the acceleration time. We define the heating time as 

tr = T /(dT/dt) . (62) 



1CT 3 10 2 10 ' 10° 10 1 10 Z 10 3 

a I A 

Fig. 3. — The comparison between heating timescale and tail- formation time-scale tp 
for different po and acceleration rates. The temperature is T = 0.016 (the corresponding 
momentum is = 0.12). The threshold value is marked by the gray horizontal line. 

The lower limit of this time can be obtained for the quasi- stationary solution for dW/ dt 
when we neglect the cooling term Sn£q in Eq. (1501) . The acceleration time characterizes a 
period required for particles to fill the non-thermal tail. Numerical calculations show that 
for £ > 2 this time is of the order of 

t F ~ a" 1 . (63) 

The quasi-stationary state (when the plasma temperature is almost constant and the 
acceleration generate prominent non-thermal "tails") can be reached only if > tp. In 
this case we can use analytical solutions presented in previous section. The ratio tj- / tp as a 
function of po is shown in Fig. [3j The threshold value of ratio tx/t-p = 1 is shown in Fig. [3] 
by the gray horizontal dashed line. The quasi- stationary state will be achieved if t T /tp is 
above the gray line. 

If the acceleration time is larger than the heating time, tr < tp, the quasi- stationary 
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state cannot be r e ached . In this case we can simplify the calculations using the trick in 



Petrosian fc East 



(120081 ) . The evolution of distribution function f{p) can be described by 



the non- stationary linear kinetic equation 



9/ + ii 2 
dt p 2 dp 



tp,N,T)f-{D (p,N,T) + D F (p)}^- 



0. 



(64) 



We can estimate the variation of temperature by the following algorithm: 

1. For a given f(t,p), estimate f(t + St,p) from Eq. ([641 ; 

2. compute N(t + St) from J °° f(t + St,p)dp; 

3. calculate W> from Eq. §U]), then W (t + St) = W Q (t) + W St, then find T(t + St) from 
Eq. (E2D; 



4. for new values of N(t + St) and T(t + St) recalculate the kinetic coefficients using 
analytical expressions for Maxwellian coefficients, see Eqs. ( 165|) & (166]) ; 

5. repeat steps 1-4. 



The analytical exp ressions for the kinetic coefficients are calculated as (see 



Petrosian fc East 



2008 



dp 
~dt 



(p,N,T) 



and references therein): 
A(p 2 + 1 



P 



erf 



4£ 

7tT 



exp 



T 



D (p,N,T) 
Here erf(z) is the error function. 



TyV + 1 /dp 

p \dt 



(p,N,T) 



(65) 



(66) 



With this method we combine the simplicity of the analytical method with the 
accuracy of the numerical method. The only problem is that this approach like any other 
semi-analytical method based on Eq. ( 15"U|) cannot be used near p = 0. 
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Now we compare the results obtained with different methods. We consider the following 
methods: 



1. Transitionless case: It is based on Eq. flM]) . The equations are integrated numerically 
using the Runge-Kutta method to obtain the evolution of the temperature T(t) and 
density N(t). This method is applicable if pi n j < po. 

2. Quasi-linear approximation: The distribution function is given as in Eq. ( 14T|) and 
( 143]) . Eqs. (|50|) and ( 152]) are used to estimate the dT/dt. The variations of the 
temperature T(i) and density N(t) are obtained using the Runge-Kutta method. This 
approximation is valid for tx > t-p- 

3. Semi-analytical method: It uses a combination of numerical solution to Eq. (I64p and 
analytical calculations of Eqs. (150]) & ( 152]) in order to estimate the variations of the 
temperature. This method can be applied to Po ^ Pt = + l) 2 — 1. 

4. Numerical method: Evolution of the distribution function is obtained by a numerical 
solution of the original non-linear equation Eq. ([7]) (for details see Appendix [B]) . This 
method is the most universal and is used to check whether the results obtained by 
methods 1-3 are correct. 



We checked our numerical program by calculations of temperature v ariations for the 



accele ration in the form Eq. ( TT5]) and for the same parameters as used by iPetrosian fc East 



(120081 ). i.e., po = 0, q = 1, £ c = 0.2 and three values of r equaled correspondingly: 2At c , 
0.18r c and 0.013r c , where t c = (Arrr^cNlnA)- 1 « 2.7 x 10 7 x (iV/10" 3 cm- 3 )- 1 yr and 
ro = e 2 /(m e c 2 ). The acceleration parameters, Dp(p), for these three cases are shown in Fig. 
[S] by the thin dashed lines. 



The result of our calculations and that of 



Petrosian fc East 



(120081 ) are shown in Fig. H 
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Fig. 4. — Comp a rison 



Petrosian fc East 



between our numerical method (solid line) and method used by 



(120081 ) (dashed line). All notations are the same as in 



Petrosian fc East 



((2008|) (see the text for details). 



by the solid and dashed lines correspondingly. One can see that despite of some discrepancy 
the results are more or less the same. 

Now we present results of calculations for the acceleration parameter Dp in the form 
( JTT1) . when p ^ 0. In all cases we take po = 0.55. Variations of po change the temperature 
variations quantitatively but not qualitatively. We choose = 2 in order to obtain the 
same momentum dependence of Dp at high energies as in Eq. ffl5|) . Note that in this case 
t-p = 2r n . Below we will use r as a characteristic timescale to compare results with those of 



Petrosian fc East 



(120081). 



For other values of ^ the results are qualitatively the same, yet lower values of <j will 
increase the amount of non-thermal particles and thus decrease the heating timescale and 
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Fig. 5. — Comparison between Dp(p) used by lPetrosian fc Eastl (120081 ) (thin dashed lines) 



Eq. (ITB]) . and in this paper (thick solid lines), Eq. (fT7|) . 
vice- versa. One can see this from Eq. (1611) . 

The calculations where performed for the three different regimes of acceleration: 

(a) heating dominates over cooling; 

(b) cooling and heating rates are of the same order; 

(c) cooling dominates over heating. 

The functions Dp(p) used for these three cases are shown in Fig. [5] by the solid lines. 

We provide calculations by the four different methods: analytical (transitionless), 
quasi-linear, semi-analytical and numerical. Temperature variations, T(t), obtained by 
these methods are shown in Figs. [fU [8] and [9] by the dashed, thin solid, thick solid and 
dotted lines, correspondingly. 
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(a) For the case of slow acceleration we take a/ A = 2.77. In this case the quasi- stationary 
approach is not valid, and only numerical and semi-analytical methods can provide 
an adequate result. Temperature variations for this case of acceleration parameter 
are shown in Fig. [61 As one can see the result of acceleration for this parameter 
a is the plasma overhea ting that is in complet e agre ement with the conclusions 



of 



Wolfe k Melia 



(J2006J) and 



Petrosian &: East 



(120081 ). The only difference is that 



th e overheating occurs for the time t w 4.5tq ~ 1.6tc which is longer than that 



of 



Petrosian &: East 



(120081 ) by a factor of 4 who obtained t ~ To. The reason is 
that for po = the acceleration generates an extended excess above the equilibrium 
Maxwellian function while for po ^ this excess is not so prominent (compare dashed 
and solid lines in Fig. |7J). Since the amount of suprathermal particles in the case 
Po = is higher than the case p > 0, it is not surprising that the plasma is overheated 
by the Coulomb losses in a shorter time when p = 0. 

(b) The case of moderate acceleration (a/ A = 11.63) is shown in Fig. [HJ One can see that 
all methods are in good agreement. At the first stage we see plasma heating, however 



the timescale is much longer than in 



Petrosian fc East 



(120081 ) . the plasma temperature 



increases by a factor of 1.3 at the moment t ~ 8 6 rn = 7tc that is almost two orders of 



magnitude higher than that of iPetrosian fc East 



(120081 ) . A prominent quasi-stationary 



power-law tail of nonthermal particles is formed by the acceleration for a much shorter 



Petrosian fc East 



(|2008|), after 



time (since tx / tj? > 1, see Fig. |3|). Moreover unlike in 
this time heating reverses to cooling. 

(c) The case of fast acceleration (a/A = 18.5) is shown in Fig. [HI The numerical method, 
the semi-analytical method and the quasi-linear approximation give almost the same 
result. For comparison we also show the calculations obtained with transitionless case 
(dashed line). We see that in spite of some difference this method provides a similar 
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Fig. 6. — The temperature evolution when heating dominates over cooling for the parameters: 
a/A = 2.77, p = 0.55, q = 2, T = 7 keV. Dotted line - numerical model, thick solid line - 
semi-analytical model, thin solid line - quasi-linear model, dashed line - transitionless model. 
Note that semi-analytical model and the numerical model almost overlap entirely. 

qualitative time variations of the temperature T. All methods demonstrate plasma 
cooling in this regime from the very beginning. The temperature of the plasma shows 
a steady decrease with time while nonthermal tails are for med rapidly (tr/ty > 1, se e 



Fig. [3]) that differs completely from the results obtained by iPetrosian & East 



fl2008f ). 



These results can be understood from Fig. [2j For the ratio a/A = 2.77 the injection 
momentum is p^ ~ 0.83, i.e., pj n j > po (similar to the upper curve of Fig. [2]). An excess 
of quasi-thermal particles is formed in the range between p and p m y Coulomb losses 
of these particles result s in effective plasma heating. As we already above-mentioned 



Petrosian &: East 



(120081 ) assumed po = that led to more extended transition region and 



more effective heating. 
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Fig. 7. — Comparison between spectra formed under influence of momentum diffusion coef- 
ficient in the form of Eq. ( fT5|) for To = 0.18tc (thin dashed lines) and those of Eq. ( fT7|) for 
a/A = 2.77, <^ = 2 and po = 0.55 (solid lines). The distribution were taken at the moment 
when the temperatures are the same. Two values of the temperature were used: T\ = 8.5 
keV, T 2 = 16.5 keV. 

In the case of a/A = 11.63, pj n j ~ 0.5 ~ po (similar to the middle curve of Fig. [2]). The 
transition region is almost negligible in this case. Therefore, plasma heating by nonthermal 
particles is insignificant which then changes into cooling. 

In the case of a/ A = 18.5, p in j ~ 0.4, i.e., p in j < p (similar to the lower curve of Fig. |2]). 
A deficit of high energy particles is formed in the thermal energy range that provides the 
effect of cooling. 

Thus, we conclude, that depending on parameters, po and a, different regimes of 
acceleration from background plasma are realized. The important inference is that stochastic 
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Fig. 8. — Same as Fig. [6] but for the case when cooling and heating processes are of the 
same order of magnitude (a/A = 11.6). 

acceleration may produce a flux of nonthermal particle without plasma overheating. 

A specific spectrum of turbulence that provides stochastic acceleration is out of the 
scope of this paper. It depends on mechanisms which excite electromagnetic fluctuations 



in an astrophysical plasma. As an example, we mention partic 



associations by a supersonic turbulence (see 



Bykov &: Toptygin 



e acceleration in 0B- 



19931 ). The momentum 



diffusion coefficient in this case has the form D(p) = D^p 2 . This acceleration is effective in 
the momentum range p > po, where the value of po is derived from r^po) — lu/c. Here tl 
is the particle Larmor radius, u is the shock velocity and I is a distance between shocks. 
Particles with p < p are not accelerated by this mechanism. 
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Fig. 9. — Same as Fig. [6]but for the case when cooling dominates over heating (a/A = 18.5). 



5. Conclusion 



We analyzed nonlinear kinetic equations describing particle stochastic (or second-order 
Fermi) acceleration from background plasma when the acceleration is non-zero for particles 
with momenta p > po- The goal of these investigations is to defi ne whether the only r esult 



of stochastic acce l eratio n is plasma overheating as concluded by 



Petrosian fc East 



Wolfe fc Melial (120061 ) and 



( 120081 ). or this acceleration can generate prominent tails of nonthermal 
particles when the plasma temperature remains almost stationary. The following results are 
obtained from our analysis: 



1. We showed that in the case of stochastic acceleration two competitive processes 
determine temperature variations of background plasma. The first one is Coulomb 
energy losses of nonthermal particles which heat the plasma. The other one is a 
run-away flux of high-energy particles from the thermal pool that leads to plasma 
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Fig. 10. — Regions of heating (below the solid line) and cooling (above the solid line) in 
the parameter space (T, a /A). The dashed lines show the evolution of systems at the same 
starting temperature but with different a /A. 

cooling. Depending on the rates of these processes the plasma may cool down or heat 
up. 

2. From numerical and analytical calculations we conclude that for a low enough 
acceleration rate the cooling process is negligible. The plasma gains much heat 
on the acceleration timescale Tq. As a result the plasma temperature rises rapidly 
while prominent non ther mal tails are not generate d, that fully confirms results of 



Wolfe fc Melial (J200J) and 



Petrosian &: East 



fl2008h . 



3. For a moderate acceleration the cooling and heating processes partly compensate each 
other. As a results the plasma temperature is quasi-stationary on a timescales much 
longer than r . In this case, the acceleration produces a nonthermal component of 
the spectrum. After a period of moderate plasma heating the process changes into 
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cooling. This reg i me do es not appear in the models of IWolfe &: Melia! ( 120061 ) and 



Petrosian fc East 



fl2008h . 



4. For a high rate of acceleration the run-away flux of thermal particles cools the plasma 
down from the very beginning. In spite of energy supply by external sources the 
plasma temperature drops down (analogue to Maxwell demon). 

5. The evolution of plasma temperature depends on the characteristic time of Coulomb 
collisions in the background plasma (the collision frequency A) and the acceleration 
frequency a. This is illustrated in Fig. [TUI where the solid line defines the border 
between heating (below the line) and cooling (above the line) regimes for quasi- 
stationary systems. It corresponds to the solution of the dT /dt = (see Eq. fl50|) ). 
Dashed lines in Fig. [10] show the evolution of plasma parameters for the same initial 
temperature but different initial value of A (a is the same for the systems). Since N 
decreases monotonically because of particle acceleration, A decreases monotonically 
accordingly (see Eq. (125]) ). 

One can see that even if the system starts from the regime of heating, sooner or 
later it changes to plasma cooling. If the evolution of the system is quasi-stationary 
the turning point of the trajectory should be located on the boundary. However we 
mention that the quasi-stationary approximation is inapplicable to low values of a /A. 



A. General Kinetic Equation 



The general equat ion for stochastic Fermi acceleration wit 



equation has the form (jLandau fc Lifshit2 



1981 



Wolfe k Melt 



i the coefficient the 



2006t ) 



dm 

dt 



_d_ 

dp 



(Al) 
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where p = v (c 2 — v 2 ) is the dimensionless particle momentum and v is the particle 
velocity. The coefficients D a » and F a are determined by Coulomb collisions of the particles, 







dp'fs 



T Z a p(p,p') 



/(p')dV 



where 



Z a p{p, p') = — ^ [w 2 5 aP - - p' a p'p + r(p a p' p + p'jpp 
r y r y id 



8vre 2 e /2 lnA 2 ^ r— — — — 

A = , r = 77— p-p/c, iy = cy)" — 1 , 7 = y 1 + p / cr 



The boundary conditions were taken in the form: a zero particle flux at p = 0: 

= 0, 



(^ + ^)^-F a /( P ) 



dp p 



and the distribution function vanishes at some p n 



p=0 



(A2) 

(A3) 
(A4) 



(A5) 



/(p. 



(A6) 



B. Numerical Method for the Nonlinear Case 



To solve Eq. (JTj) numerically we use the Crank-Nicolson finite difference method. To 
estimate the kinetic coefficients flA2j) we use Simpson's integration rule, 



D = Zf. 



Z'f. 



(Bl) 



where vectors f, D and F are corresponding discrete versions of f(p), D c (p, f) and 
(dp/dt) c (p, /). Matrices Z and Z' are obtained by applying Simpson's rule to Eq. ( 1A2j) . 
The discrete vers i on of Eq. (EJ) at t = (t n + t n+ i)/2 and p = pj looks like (see e.g., 



Park fc Petrosian 



(119961 ) and references therein) 

fn+l,j fn,j 



1 S n+\,j+\ S n+\,j- 



At 



Pi 



(B2) 
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where At = t n+ i — t n and Apj = (pj+i — pj_i)/2 are steps of the grid and the flux is 
expressed according to the Crank-Nicolson rule: 

/n+lj+l — fn+l,j 



I„ 2 
20 



i+i 



Ap 



D 



fn,j+l fn, 



n,j+{ 



— F, 



n, j+ lfn, j+ l 



+ 



(B3) 



D 



fn 



2{fn,j + fn,j+l) i 
Pj+1 - Pj ■ 



(B4) 
(B5) 
(B6) 
(B7) 



The boundary condition at p = imply that S n+ i_i = 0. The boundary condition at 

V = Pmax IS fiPmax) = 0. 



After discretization we arrive at the non-linear system of equations, 



fn+l - fn = -4(fn+l)fn+l + A%)f n 



(B8) 



where A(f) is a tridiagonal matrix corresponding to the differential operator in RHS of 
Eq. (1B2I) . According to Eq. (1B1I) . A(f) is a linear function of f and Eq. (1B8j) is a system of 
quadratic equations. 

To avoid calculations of Jacobian matrix we do not apply Newton's method and 
utilize a simple iteration method instead. However the iteration method based on Eq. (IB8I) 
converges very slowly. We rewrite the iteration step in the following form 



(E - A(f* +1 )) £g = A(f a % + f n 



(B9) 



where k is the number of iteration and E is a unit matrix. The system of linear equations 
is solved using tridiagonal matrix algorithm. Iteration in the form Eq. (1B9[) shows fast 



convergence if the temperature of the Maxwellian distribution does not change significantly 
between t n and t n+1 . 

Since the Crank-Nicolson method may be affected by numerical oscillations we also 
use l ess precise and more rob ust backward Euler method (simple fully implicit method 



from 



Park &: Petrosian 



19961 ) . The backward Euler method turns out to be useful for 



non-thermal tails of low magnitude when the value of a is low and the value of po is high. 

The discretization in momentum space is tricky since we need to provide a good 
resolution for Maxwellian distribution and transitional region as well as calculate the 
nonthermal tail at high energies. We tried two possible ways to reduce the number of grid 
points. The first one is to split the momentum axis into sub-domains and join them using 
continuity of the distribution function and particle flux. The second way is to use the 
logarithmic grid by introducing new variable q = log(p). Both methods gives almost the 
same results. 
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